// The limits of ‘Western’ supply chain sustainability governance to halt deforestation
// Lead Author: Christoph Kubitza
//Email: Christoph.Kubitza@giga-hamburg.de
// Date: 12.4.2025


from qgis.core import (
    QgsProject,
    QgsVectorLayer,
    QgsProcessingFeedback,
    QgsApplication
)
from qgis.analysis import QgsNativeAlgorithms
import processing

##Generate grid cells

# Initialize QGIS Application (if running as standalone script)
qgs = QgsApplication([], False)
qgs.initQgis()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Define file paths
base_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/Base Layer/base_layer500m.shp"
concession_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Concession Maps/oilpalm_deals_merged_55.shp"
output_intersection_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final"

# Load input layers
base_layer = QgsVectorLayer(base_layer_path, "Base Layer", "ogr")
concession_layer = QgsVectorLayer(concession_layer_path, "Concession Layer", "ogr")

# Validate layers
if not base_layer.isValid():
    raise Exception("Base layer failed to load.")
if not concession_layer.isValid():
    raise Exception("Concession layer failed to load.")

QgsProject.instance().addMapLayer(base_layer)
QgsProject.instance().addMapLayer(concession_layer)

# Define intersection parameters
params = {
    'INPUT': concession_layer_path,  # Path to the concession layer
    'INPUT_FIELDS': [],  # No specific fields from INPUT
    'OVERLAY': base_layer_path,  # Path to the base layer
    'OVERLAY_FIELDS': [],  # No specific fields from OVERLAY
    'OVERLAY_FIELDS_PREFIX': '',  # No prefix for OVERLAY fields
    'OUTPUT': 'C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final'  # Path to the output file
}

# Run the intersection process
result = processing.run("native:intersection", params)

# Load the resulting intersection layer
intersection_layer = QgsVectorLayer(result['OUTPUT'], "Intersection Layer", "ogr")

# Validate the resulting layer
if not intersection_layer.isValid():
    raise Exception("Intersection layer failed to load.")
else:
    QgsProject.instance().addMapLayer(intersection_layer)
    print("Intersection layer successfully added to the project.")

# Define the output file path
output_file_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C.shp"

# Save the intersection_layer as a shapefile
error = QgsVectorFileWriter.writeAsVectorFormat(
    intersection_layer,
    output_file_path,
    "UTF-8",  # Encoding
    intersection_layer.crs(),  # Coordinate Reference System
    "ESRI Shapefile"  # Format
)

if error == QgsVectorFileWriter.NoError:
    print(f"Intersection layer successfully saved to {output_file_path}")
else:
    print(f"Error saving layer: {error}")




#Split grid cells that not are contiguous (otherwise cannot use centroids later on)
from qgis.core import (
    QgsApplication,
    QgsProcessingFeedback
)
from qgis.analysis import QgsNativeAlgorithms
from qgis import processing
from processing.core.Processing import Processing
import os

# Initialize Processing and add provider
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Define file paths
input_layer_path = r"C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C.shp"
output_layer_path = r"C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp.shp"

# Check if input file exists before proceeding
if not os.path.exists(input_layer_path):
    print(f"Error: Input file does not exist: {input_layer_path}")
    qgs.exitQgis()
    exit()

# Run "Multipart to Singleparts"
params = {
    'INPUT': input_layer_path,
    'OUTPUT': output_layer_path
}

try:
    result = processing.run("native:multiparttosingleparts", params, feedback=QgsProcessingFeedback())
    if result and os.path.exists(output_layer_path):
        print(f"Multipart to Singleparts conversion successful! Output saved to: {output_layer_path}")
    else:
        print("Error: Processing failed.")
except Exception as e:
    print(f"Processing error: {e}")




## Intersect layer with peatlands and create unique IDs
import processing

# Define file paths
input_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp.shp"
peatland_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Peatland/Indonesia_peat_lands.shp"
output_intersection_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp_union.shp"

# Load the layers
input_layer = QgsVectorLayer(input_layer_path, "Input Layer", "ogr")
peatland_layer = QgsVectorLayer(peatland_layer_path, "Peatland Layer", "ogr")

# Check if layers are valid
if not input_layer.isValid() or not peatland_layer.isValid():
    raise ValueError("Error loading one or both layers")

print("Both layers loaded successfully")

# Ensure both layers have the same CRS
if input_layer.crs() != peatland_layer.crs():
    print("Reprojecting input layer to match peatland layer CRS...")
    input_layer = processing.run("native:reprojectlayer", {
        "INPUT": input_layer,
        "TARGET_CRS": peatland_layer.crs(),
        "OUTPUT": "memory:"
    })["OUTPUT"]

# Step 1: Get intersection (splitted part)
print("Running intersection...")
intersection_result = processing.run("native:intersection", {
    "INPUT": input_layer,
    "OVERLAY": peatland_layer,
    "OUTPUT": "memory:"
})["OUTPUT"]

# Step 2: Get difference (remaining part of polygons)
print("Running difference to extract non-intersecting parts...")
difference_result = processing.run("native:difference", {
    "INPUT": input_layer,
    "OVERLAY": peatland_layer,
    "OUTPUT": "memory:"
})["OUTPUT"]

# Step 3: Merge the split and non-intersecting parts back together
print("Merging results...")
final_result = processing.run("native:mergevectorlayers", {
    "LAYERS": [intersection_result, difference_result],
    "OUTPUT": output_intersection_path
})["OUTPUT"]

if final_result:
    print(f"Processing completed successfully! Output saved at: {output_intersection_path}")
else:
    print("Processing failed.")


output_inters_layer = QgsVectorLayer(output_intersection_path, "Output Intersect Layer", "ogr")
unique_id_field = QgsField("unique_id", QVariant.Int)
output_inters_layer.dataProvider().addAttributes([unique_id_field])
output_inters_layer.updateFields()

expression = QgsExpression('@row_number')
context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(output_inters_layer))

with edit(output_inters_layer):
    for i, feature in enumerate(output_inters_layer.getFeatures()):
        context.setFeature(feature)
        feature["unique_id"] = i + 1
        output_inters_layer.updateFeature(feature)

# Save changes
output_inters_layer.commitChanges()

print("Process completed. Intersection layer with unique IDs is ready.")



##AREA
# Import required modules
from qgis.core import (
    QgsVectorLayer,
    QgsField,
    QgsProject,
    QgsCoordinateReferenceSystem,
    QgsCoordinateTransform,
    QgsVectorFileWriter
)
from PyQt5.QtCore import QVariant

# Define input and output file paths
input_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp_union.shp"
output_reprojected_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp_union_proj.shp"

# Load the original layer
input_layer = QgsVectorLayer(input_path, "Original Layer", "ogr")

# Check if layer is valid
if not input_layer.isValid():
    print("Layer failed to load!")
else:
    print("Layer loaded successfully!")

# Define the target CRS (EPSG:32752 - WGS 84 / UTM zone 52S)
target_crs = QgsCoordinateReferenceSystem("EPSG:32752")

#Here we did manually reproject in case it did not work
# Reproject the layer
reprojected_layer = QgsVectorFileWriter.writeAsVectorFormat(
    input_layer, output_reprojected_path, "UTF-8", target_crs, "ESRI Shapefile"
)

if reprojected_layer == QgsVectorFileWriter.NoError:
    print(f"Reprojection successful! Saved to {output_reprojected_path}")
else:
    print("Reprojection failed!")

# Load the reprojected layer
layer = QgsVectorLayer(output_reprojected_path, "Reprojected Layer", "ogr")

# Start editing
layer.startEditing()

# Check if the area field already exists
field_names = [field.name() for field in layer.fields()]
if "area_ha" not in field_names:
    layer.dataProvider().addAttributes([QgsField("area_ha", QVariant.Double)])
    layer.updateFields()

# Get the index of the new field
area_idx = layer.fields().indexFromName("area_ha")

# Iterate over features and calculate the area
for feature in layer.getFeatures():
    geom = feature.geometry()
    if geom and geom.isGeosValid():
        area_m2 = geom.area()  # Area in square meters (since CRS is in meters)
        area_m2 = area_m2   
        feature.setAttribute(area_idx, area_m2)
        layer.updateFeature(feature)

# Save the changes
layer.commitChanges()
print("Area calculation complete!")

# Optionally add the layer to QGIS
QgsProject.instance().addMapLayer(layer)


#check if everything went well using this command in the field calculator: area(transform($geometry, @layer_crs, 'EPSG:32752'))



## INTERSECTION WITH RASTER FILES ZONAL STAT
from qgis.core import (
    QgsVectorLayer,
    QgsProject,
    QgsProcessingFeatureSourceDefinition
)
from qgis.analysis import QgsNativeAlgorithms  # Ensure native QGIS algorithms are available
import processing
from qgis.analysis import QgsZonalStatistics


# Define file paths
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_sp_union_proj.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats.shp"

# Load intersection layer
intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

# GAEZ Load raster layer
rasterLayergaez = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Oilpalm Yield GAEZ/OP national medium input att yield/attyield_m.tif')
zoneStat = QgsZonalStatistics(intersection_layer, rasterLayergaez, 'gaez_', 1)
zoneStat.calculateStatistics(None)
# Initialize QGIS processing algorithms (if not already done)
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Acessability Load raster layer
rasterLayeraccess = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Accessability/cost_distance_IND_roads.tif')
zoneStat = QgsZonalStatistics(intersection_layer, rasterLayeraccess, 'road_', 1)
zoneStat.calculateStatistics(None)
# Initialize QGIS processing algorithms (if not already done)
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# SRTM Load raster layer
rasterLayersrtm = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Altitude SRTM/IND_srtm_250m_wgs84.tif')
zoneStat = QgsZonalStatistics(intersection_layer, rasterLayersrtm, 'srtm_', 1)
zoneStat.calculateStatistics(None)
# Initialize QGIS processing algorithms (if not already done)
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Forest cover Load raster layer
rasterLayerhansen = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Forest cover Hansen/Raster layer/CanopyCover10_2000.tif')
zoneStat = QgsZonalStatistics(intersection_layer, rasterLayerhansen, 'ccover_', 1)
zoneStat.calculateStatistics(None)
# Initialize QGIS processing algorithms (if not already done)
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

error = QgsVectorFileWriter.writeAsVectorFormat(
    intersection_layer,
    output_layer_path,
    "UTF-8",  # Encoding
    intersection_layer.crs(),  # Coordinate Reference System
    "ESRI Shapefile"  # Format
)




## INTERSECTION WITH RASTER FILES ZONAL HIST
from qgis.core import QgsApplication, QgsRasterLayer, QgsVectorLayer
from qgis.analysis import QgsNativeAlgorithms
from qgis import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Paths to input files
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo.shp"


# Ensure layers are loaded
# Load intersection layer
intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

#Fire data
rasterLayer = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Forest fire/SEA-AUS_fire_forest_loss_2001-22_annual.tif')

hist_result = processing.run("native:zonalhistogram", {
    'INPUT_RASTER': rasterLayer,
    'RASTER_BAND': 1,
    'INPUT_VECTOR': intersection_layer,
    'COLUMN_PREFIX': 'fire_',
    'OUTPUT': output_layer_path
})

#Danylo data
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo2.shp"
rasterLayer = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND OilPalm Danylo/oilpalm8417.tif')

intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

hist_result = processing.run("native:zonalhistogram", {
    'INPUT_RASTER': rasterLayer,
    'RASTER_BAND': 1,
    'INPUT_VECTOR': intersection_layer,
    'COLUMN_PREFIX': 'op_year',
    'OUTPUT': output_layer_path
})

#Gaveau data
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo2.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo3.shp"
rasterLayer = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND OilPalm Gaveau/Industrial_Oilpalm_Kalimantan_Papua_30m_wgs84.tif')

intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

hist_result = processing.run("native:zonalhistogram", {
    'INPUT_RASTER': rasterLayer,
    'RASTER_BAND': 1,
    'INPUT_VECTOR': intersection_layer,
    'COLUMN_PREFIX': 'gav_',
    'OUTPUT': output_layer_path
})

#JRC data
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo3.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo4.shp"
rasterLayer = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Forest cover JRC/DefYear2023_JRC.tif')

intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

hist_result = processing.run("native:zonalhistogram", {
    'INPUT_RASTER': rasterLayer,
    'RASTER_BAND': 1,
    'INPUT_VECTOR': intersection_layer,
    'COLUMN_PREFIX': 'jrc_',
    'OUTPUT': output_layer_path
})

#Hansen data
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo4.shp"
output_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5.shp"
rasterLayer = QgsRasterLayer('C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Forest cover Hansen/Raster layer/treeLoss0021.tif')

intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    raise Exception("Failed to load intersection layer.")
QgsProject.instance().addMapLayer(intersection_layer)

hist_result = processing.run("native:zonalhistogram", {
    'INPUT_RASTER': rasterLayer,
    'RASTER_BAND': 1,
    'INPUT_VECTOR': intersection_layer,
    'COLUMN_PREFIX': 'treeloss',
    'OUTPUT': output_layer_path
})





## PEATLAND LAYER
from qgis.core import QgsApplication, QgsRasterLayer, QgsVectorLayer, QgsProject
from qgis.analysis import QgsNativeAlgorithms
from qgis import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Define file paths
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5.shp"
peatland_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Peatland/Projected/Indonesia_peat_lands_proj.shp"
output_projected_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84.shp"
output_nn_join_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_nn_join.shp"
final_output_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_peat.shp"

# Load layers
intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
peatland_layer = QgsVectorLayer(peatland_layer_path, "Peatland Layer", "ogr")

if not intersection_layer.isValid() or not peatland_layer.isValid():
    raise Exception("One or more layers failed to load.")

QgsProject.instance().addMapLayer(intersection_layer)
QgsProject.instance().addMapLayer(peatland_layer)

# Reproject the intersection layer to EPSG:32752 (WGS 84 / UTM zone 52S)
reproject_params = {
    'INPUT': intersection_layer,
    'TARGET_CRS': 'EPSG:32752',
    'OUTPUT': output_projected_layer_path
}
reproject_result = processing.run("native:reprojectlayer", reproject_params)

projected_layer = QgsVectorLayer(output_projected_layer_path, "Projected Intersection Layer", "ogr")
if not projected_layer.isValid():
    raise Exception("Reprojected layer failed to load.")
QgsProject.instance().addMapLayer(projected_layer)



## Create Centroids
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84.shp"
output_centroid_layer = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids.shp"

# Load the layer
intersection_layer = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
if not intersection_layer.isValid():
    print("Failed to load the intersection layer!")
else:
    # Create a new layer for centroids
    centroid_layer = QgsVectorLayer("Point?crs=" + intersection_layer.crs().authid(), "Centroids", "memory")
    centroid_layer_provider = centroid_layer.dataProvider()
    
    # Copy fields from the original layer
    centroid_layer_provider.addAttributes(intersection_layer.fields())
    centroid_layer.updateFields()

    # Generate centroids
    features = []
    for feature in intersection_layer.getFeatures():
        centroid = feature.geometry().centroid()
        new_feature = QgsFeature()
        new_feature.setGeometry(centroid)
        new_feature.setAttributes(feature.attributes())  # Keep original attributes
        features.append(new_feature)

    # Add centroids to the new layer
    centroid_layer_provider.addFeatures(features)
    centroid_layer.updateExtents()

    # Save to shapefile
    QgsVectorFileWriter.writeAsVectorFormat(centroid_layer, output_centroid_layer, "UTF-8", centroid_layer.crs(), "ESRI Shapefile")
    print(f"Centroid layer saved at: {output_centroid_layer}")

# Add the layer to the QGIS interface
QgsProject.instance().addMapLayer(QgsVectorLayer(output_centroid_layer, "Centroids", "ogr"))




##Measure nearest distance 
#Our peatlayer is already in EPSG:32752 - WGS 84
from qgis.core import (
    QgsProject, QgsVectorLayer, QgsFeature, QgsGeometry,
    QgsPointXY, QgsVectorDataProvider, QgsField
)
from PyQt5.QtCore import QVariant

# Define paths
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids.shp"
peatland_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Peatland/Projected/Indonesia_peat_lands_proj.shp"

# Load layers
p_lyr = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
l_lyr = QgsVectorLayer(peatland_layer_path, "Peatland Layer", "ogr")

# Check if layers are valid
if not p_lyr.isValid() or not l_lyr.isValid():
    print("Error loading one or both layers")
else:
    QgsProject.instance().addMapLayer(p_lyr)
    QgsProject.instance().addMapLayer(l_lyr)

# Create memory layer for distances
d_lyr = QgsVectorLayer('LineString?field=distance:float', 'dist', 'memory')
QgsProject.instance().addMapLayer(d_lyr)

# Get data provider
prov = d_lyr.dataProvider()

feats = []
distances = {}

for p in p_lyr.getFeatures():
    p_geom = p.geometry()
    
    if not p_geom.isEmpty() and not p_geom.isMultipart():
        p_point = p_geom.asPoint()
        minDist = float("inf")
        minDistPoint = None

        # Find closest point on the peatland layer
        for l in l_lyr.getFeatures():
            l_geom = l.geometry()
            if not l_geom.isEmpty():
                dist = p_geom.distance(l_geom)  # True minimum distance

                if dist < minDist:
                    minDist = dist
                    minDistPoint = l_geom.nearestPoint(p_geom).asPoint()

        if minDistPoint:
            feat = QgsFeature()
            feat.setGeometry(QgsGeometry.fromPolyline([
                QgsPoint(p_point), QgsPoint(minDistPoint.x(), minDistPoint.y())
            ]))
            feat.setAttributes([minDist])
            feats.append(feat)
            distances[p.id()] = minDist  # Store correct distance

# Add features to the layer
prov.addFeatures(feats)
d_lyr.updateExtents()
d_lyr.triggerRepaint()

# Add distance field to p_lyr if not present
p_lyr_provider = p_lyr.dataProvider()
if "distance" not in [field.name() for field in p_lyr.fields()]:
    p_lyr_provider.addAttributes([QgsField("distance", QVariant.Double)])
    p_lyr.updateFields()

# Update p_lyr with calculated distances
p_lyr.startEditing()
for p in p_lyr.getFeatures():
    if p.id() in distances:
        p.setAttribute(p_lyr.fields().indexFromName("distance"), distances[p.id()])
        p_lyr.updateFeature(p)
p_lyr.commitChanges()

output_peat_layer = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids_p.shp"

QgsVectorFileWriter.writeAsVectorFormat(p_lyr, output_peat_layer, "UTF-8", p_lyr.crs(), "ESRI Shapefile")
print(f"Centroid layer saved at: {output_peat_layer}")



#Measure nearest distance from inside the peatlands
#Our peatlayer is already in EPSG:32752 - WGS 84
from qgis.core import (
    QgsProject, QgsVectorLayer, QgsFeature, QgsGeometry,
    QgsPointXY, QgsVectorDataProvider, QgsField, QgsVectorFileWriter
)
from PyQt5.QtCore import QVariant

# Define paths
intersection_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids_p.shp"
peatland_layer_path = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND Peatland/Projected/Indonesia_peat_lands_outside_proj.shp"

# Load layers
p_lyr = QgsVectorLayer(intersection_layer_path, "Intersection Layer", "ogr")
l_lyr = QgsVectorLayer(peatland_layer_path, "Peatland Layer", "ogr")

# Check if layers are valid
if not p_lyr.isValid() or not l_lyr.isValid():
    print("Error loading one or both layers")
else:
    QgsProject.instance().addMapLayer(p_lyr)
    QgsProject.instance().addMapLayer(l_lyr)

# Create memory layer for distances
d_lyr = QgsVectorLayer('LineString?field=dist_in:float', 'dist', 'memory')
QgsProject.instance().addMapLayer(d_lyr)

# Get data provider
prov = d_lyr.dataProvider()

feats = []
distances = {}

for p in p_lyr.getFeatures():
    p_geom = p.geometry()
    
    if not p_geom.isEmpty() and not p_geom.isMultipart():
        p_point = p_geom.asPoint()
        minDist = float("inf")
        minDistPoint = None

        # Find closest point on the peatland layer
        for l in l_lyr.getFeatures():
            l_geom = l.geometry()
            if not l_geom.isEmpty():
                dist = p_geom.distance(l_geom)  # True minimum distance

                if dist < minDist:
                    minDist = dist
                    minDistPoint = l_geom.nearestPoint(p_geom).asPoint()

        if minDistPoint:
            feat = QgsFeature()
            feat.setGeometry(QgsGeometry.fromPolyline([
                QgsPoint(p_point), QgsPoint(minDistPoint.x(), minDistPoint.y())
            ]))
            feat.setAttributes([minDist])
            feats.append(feat)
            distances[p.id()] = minDist  # Store correct distance

# Add features to the layer
prov.addFeatures(feats)
d_lyr.updateExtents()
d_lyr.triggerRepaint()

# Add dist_in field to p_lyr if not present
p_lyr_provider = p_lyr.dataProvider()
if "dist_in" not in [field.name() for field in p_lyr.fields()]:
    p_lyr_provider.addAttributes([QgsField("dist_in", QVariant.Double)])
    p_lyr.updateFields()

# Update p_lyr with calculated distances
p_lyr.startEditing()
for p in p_lyr.getFeatures():
    if p.id() in distances:
        p.setAttribute(p_lyr.fields().indexFromName("dist_in"), distances[p.id()])
        p_lyr.updateFeature(p)
p_lyr.commitChanges()

output_peat_layer = "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids_p2.shp"

QgsVectorFileWriter.writeAsVectorFormat(p_lyr, output_peat_layer, "UTF-8", p_lyr.crs(), "ESRI Shapefile")
print(f"Centroid layer saved at: {output_peat_layer}")



##Merge RSPO layer
from qgis.core import (
    QgsApplication,
    QgsVectorLayer,
    QgsProcessingFeedback
)
from qgis.analysis import QgsNativeAlgorithms
from qgis import processing
from processing.core.Processing import Processing
import os

# Initialize QGIS application (for standalone scripts)
QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis", True)  # Update if needed
qgs = QgsApplication([], False)
qgs.initQgis()

# Initialize Processing and add provider
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

# Define file paths
intersection_layer_path = r"C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/01_raw/IND RSPO/oilpalm_deals_merged_53_RSPO.shp"
input_layer_path = r"C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids_p2.shp"
output_layer_path = r"C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_C_stats_histo5_wgs84_centroids_p2_r.shp"

# Check if files exist before proceeding
if not os.path.exists(intersection_layer_path) or not os.path.exists(input_layer_path):
    print("Error: One or more input files do not exist.")
    qgs.exitQgis()
    exit()

# Load input layers
input_layer = QgsVectorLayer(input_layer_path, "input_layer", "ogr")
intersection_layer = QgsVectorLayer(intersection_layer_path, "intersection_layer", "ogr")

# Ensure layers are loaded correctly
if not input_layer.isValid() or not intersection_layer.isValid():
    print("Error: One or both layers failed to load.")
    qgs.exitQgis()
    exit()
else:
    print("Layers loaded successfully.")

# Run "Join attributes by location" with a processing feedback object
params = {
    'INPUT': input_layer_path,  # Use file path instead of layer object
    'JOIN': intersection_layer_path,
    'PREDICATE': [0],  # Intersects
    'JOIN_FIELDS': ['rspo_deal'],  # Only merge this field
    'METHOD': 1,  # Keep all records
    'DISCARD_NONMATCHING': False,  # Keep all input features
    'PREFIX': '',  # No prefix for joined fields
    'OUTPUT': output_layer_path
}

try:
    result = processing.run("native:joinattributesbylocation", params, feedback=QgsProcessingFeedback())
    if result and os.path.exists(output_layer_path):
        print(f"Output saved to: {output_layer_path}")
    else:
        print("Error: Processing failed.")
except Exception as e:
    print(f"Processing error: {e}")




##EXPLANATION OF PYTHON CODE

Simplified Workflow in QGIS v2

1. Open the layer of 500m pixels (C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\Base Layer\base_layer500m.shp)

2. Generate a intersection of C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\Base\base_layer500m.shp with C:\Users\christoph.kubitza\OneDrive - GIGA\Projects\CK Oil palm investors\Data\IND Concession Maps/oilplam_deals_merged_55.shp. 
Split polygons along peatland layer

3. Generate unique ID for pixels with @row_number based on the explantion on https://www.northrivergeographic.com/qgis-unique-field-values/

4. Calculate area of each pixel with @area

5. Save file in "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_Concession.shp

5. Use file in "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_Concession.shp

6. Use QGIS zonal statistics command based on layer "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_Concession.shp using the following files:
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Oilpalm Yield GAEZ\OP national medium input att yield\attyield_m.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Accessability\cost_distance_IND_roads.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Altitude SRTM\IND_srtm_250m_wgs84.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Forest cover Hansen\Raster layer\forestcover2000.tif

7. Use QGIS zonal histogram command based on layer  "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_Bl_Concession.shp using the following files:
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Forest cover Hansen\Raster layer\treeLoss0021.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Forest cover JRC\DefYear2023_JRC.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Forest fire\SEA-AUS_fire_forest_loss_2001-22_annual.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND OilPalm Danylo\oilpalm8417.tif
C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND OilPalm Gaveau\Industrial_Oilpalm_Kalimantan_Papua_30m_wgs84.tif

8. Save file in "C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_rasters.shp

9. Project to CRS EPSG:32752 - WGS 84 / UTM zone 52S - Projected and use QGIS command NN Join (Centroids of C:/Users/christoph.kubitza/OneDrive - GIGA/Working Paper/CK Oil palm investors/Data/03_temp/final/Intersection_rasters.shp) for the peatland layer C:\Users\christoph.kubitza\OneDrive - GIGA\Working Paper\CK Oil palm investors\Data\01_raw\IND Peatland\Projected\Indonesia_peat_lands_proj.shp

10. Same process for the outside peatlands

11. Attributes join with RSPO layer




